#ifndef DDMCMC_Vector
#define DDMCMC_Vector

#include "Matrix.h"


template<class Tp>
class mcmcVector : public mcmcMatrix<Tp>{
public:
	// construction & destruction
	mcmcVector(){};
	mcmcVector(const mcmcVector<Tp>& a);
	mcmcVector(const int row, const int col);
	mcmcVector(const int size){ SetRowVector(size); };
	~mcmcVector(){};

	// operators
	mcmcVector<Tp>&		operator=(const mcmcVector<Tp>& x);
	__forceinline Tp&	operator()(const int index);
	__forceinline Tp&	operator()(const int row, const int col);

	__forceinline mcmcVector operator+(const mcmcVector<Tp> &x);
	__forceinline mcmcVector operator-(const mcmcVector<Tp> &x);
	__forceinline mcmcVector operator*(const mcmcMatrix<Tp> &x);
	__forceinline mcmcVector operator*(const mcmcVector<Tp> &x);

	friend std::ostream &operator<< (std::ostream &strm, mcmcVector& vec){ return strm << vec.data << endl; };
	void operator>>(std::ofstream &strm);
	void operator<<(std::ifstream &strm);

	__forceinline mcmcVector	operator+(const Tp y);
	__forceinline mcmcVector	operator-(const Tp y);
	__forceinline mcmcVector	operator*(const Tp y);
	__forceinline mcmcVector	operator/(const Tp y);

	// operations
	void	SetRowVector(const int size);
	void	SetColVector(const int size);
	bool	RowVector(void);
	bool	ColVector(void);
	bool	ExpandTo(const int size);


	double	Norm(const double l);
	Tp		Projection(const mcmcVector &a);
	Tp		DotProduct(const mcmcVector &a);

	void	CopyContent(const mcmcVector<Tp> &b, const int num);
	virtual mcmcVector  Conv1D(const mcmcVector& kern, const bool bcircular);	// circular convolution in 1D
};


////////////////////////////////////////////
// construction & destruction
template<class Tp>
mcmcVector<Tp>::mcmcVector(const mcmcVector<Tp>& a){
	(*this) = a;
}

template<class Tp>
mcmcVector<Tp>::mcmcVector(const int row, const int col){
	SetDimension(row, col);
}

////////////////////////////////////////////
// operators

template<class Tp>
mcmcVector<Tp>& mcmcVector<Tp>::operator=(const mcmcVector<Tp>& x){
	//beware of self assignments : x = x
	if (this != &x){
		if (rows() != x.rows() || cols() != x.cols()){
			SetDimension(x.rows(), x.cols());
		}

		data = x.data.clone();
	}
	return *this;
}

template<class Tp>
void mcmcVector<Tp>::operator>>(ofstream &strm)
{
	// dimensions
	strm.write((char*)&nrow, sizeof(nrow));
	strm.write((char*)&ncol, sizeof(ncol));

	// data
	strm.write((char*)data.data, sizeof(Tp)*size());
}

template<class Tp>
void mcmcVector<Tp>::operator<<(ifstream &strm)
{
	int row_read, col_read;

	MatrixFree();

	// dimensions
	strm.read((char*)&row_read, sizeof(nrow));
	strm.read((char*)&col_read, sizeof(ncol));

	if (row_read>0 && col_read>0)
		SetDimension(row_read, col_read);
	// data
	strm.read((char*)data.data, sizeof(Tp)*size());
}

template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator+ (const mcmcVector<Tp>& x){
	mcmcVector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary matrix
	temp.data = data + x.data;

	return temp; //return the new matrix
}

template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator+ (const Tp y){
	mcmcVector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary matrix
	temp.data = data + y;

	return temp; //return the matrix
}

template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator - (const mcmcVector<Tp>& x){
	mcmcVector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary mcmcVector
	temp.data = data - x.data;

	return temp; //return the new mcmcVector
}

template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator - (const Tp y){
	mcmcVector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary mcmcVector
	temp.data = data - y;

	return temp; //return the mcmcVector
}

template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator * (const mcmcVector<Tp>& x){
	mcmcVector<Tp>	temp(nrow, x.ncol); //intialize the temporary mcmcVector
	cv::Mat a, b, c;
	if (cv::DataType<Tp>::type != CV_64F){
		data.convertTo(a, CV_64F);
		x.data.convertTo(b, CV_64F);
		c = a*b;
		c.convertTo(c, cv::DataType<Tp>::type);
		temp.data = c;
		return temp;
	}
	else{
		a = data.clone();
		b = x.data.clone();
		c = a*b;
		temp.data = c;
		return temp;
	}
}

template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator * (const mcmcMatrix<Tp>& x){
	Tp			sum;
	mcmcVector<Tp>	temp;

	temp.SetDimension(rows(), x.cols()); //intialize the temporary Vector
	for (int i = 0; i < nrow; i++)
	for (int j = 0; j < x.cols(); j++)
	{
		sum = 0;
		for (int k = 0; k < ncol; k++)
			sum = sum + Data(i, k)*x.Data(k, j); //multiply

		temp(i, j) = sum;
	}

	return temp; //return the result of multiplication
}

template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator * (const Tp y){
	mcmcVector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary mcmcVector
	temp.data = data*y;

	return temp; //return the mcmcVector
}


template<class Tp>
__forceinline mcmcVector<Tp> mcmcVector<Tp>::operator / (const Tp y)
{  //multiply by a scalar constant
	mcmcVector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary mcmcVector
	temp.data = data/y;

	return temp; //return the mcmcVector
}

template<class Tp>
__forceinline Tp& mcmcVector<Tp>::operator ()(const int row, const int col){
	return data.at<Tp>(row, col);
}

template<class Tp>
__forceinline Tp& mcmcVector<Tp>::operator ()(const int index){
	return data.at<Tp>(index / ncol, index % ncol);
}

template<class Tp>
void mcmcVector<Tp>::SetRowVector(const int size)
{
	SetDimension(1, size);
}

template<class Tp>
void mcmcVector<Tp>::SetColVector(const int size)
{
	SetDimension(size, 1);
}

template<class Tp>
bool mcmcVector<Tp>::ExpandTo(const int size)
{
	if (RowVector())
		return mcmcMatrix<Tp>::ExpandTo(1, size);
	else
		return mcmcMatrix<Tp>::ExpandTo(size, 1);
}

template<class Tp>
bool mcmcVector<Tp>::RowVector(void)
{
	if (cols()>0 && rows() == 1)
		return true;
	else
		return false;
}

template<class Tp>
bool mcmcVector<Tp>::ColVector(void)
{
	if (cols() == 1 && rows()>0)
		return true;
	else
		return false;
}

template<class Tp>
Tp mcmcVector<Tp>::Projection(const mcmcVector<Tp> &a){
	Tp sum = 0;
	for (int i = 0; i<size(); i++)
		sum = sum + Data(i)*a.Data(i);
	return sum;
}

template<class Tp>
Tp mcmcVector<Tp>::DotProduct(const mcmcVector<Tp> &a){
	Tp sum = 0;
	for (int i = 0; i<size(); i++)
		sum = sum + Data(i)*a.Data(i);
	return sum;
}

template<class Tp>
double mcmcVector<Tp>::Norm(const double l)
{
	int		i;
	double	sum = 0.0;

#ifdef _DEBUG
	assert(cols()>0);
#endif

	for (i = 0; i<rows(); i++)
		sum = sum + pow(fabs((double)Data(i)), l);
	sum = pow(sum, 1.0 / l);
	return sum;
	/* Jiadi
	if (l == 1){
		return cv::norm(data, cv::NORM_L1);
	}else if (l == 2){
		return cv::norm(data, cv::NORM_L2);
	}else{
		int		i;
		double	sum = 0.0;

		for (i = 0; i<size(); i++)
			sum = sum + pow(fabs((double)Data(i)), l);
		sum = pow(sum, 1.0 / l);
		return sum;
	}
	*/
}

template<class Tp>
void mcmcVector<Tp>::CopyContent(const mcmcVector<Tp> &b, const int num)
{
	for (int i = 0; i<num; i++)
		(*this)(i) = b.Data(i);
}

// circular convolution in 1D
template<class Tp>
mcmcVector<Tp> mcmcVector<Tp>::Conv1D(const mcmcVector& kern, const bool bcircular)
{
	mcmcVector<Tp>	I;
	int		i, k, m, nsize, nsize_kernel;

	I.SetDimension(rows(), cols());
	I.Zero();

	nsize = size();
	nsize_kernel = kern.size();

	if (bcircular)
	{
		// circular convolution
		for (i = 0; i<nsize; i++)
		{
			for (k = 0; k<nsize_kernel; k++)
			{
				m = (i - (k - (nsize_kernel - 1) / 2) + nsize * 100) % nsize;

				I(i) = I(i) + Data(m)*kern.Data(k);
			}
		}
	}
	else
	{
		for (i = 0; i<nsize; i++)
		{
			for (k = 0; k<nsize_kernel; k++)
			{
				m = i - (k - (nsize_kernel - 1) / 2);

				if (m >= 0 && m<nsize)
					I(i) = I(i) + Data(m)*kern.Data(k);
			}
		}
	}

	return I;
}


#endif
